gusucode.com > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM源码程序 > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM\stprtool\quadrat\quad2d.m
function [X1,Y1,X2,Y2]=quad2d(A,B,C,win,interp) % QUAD2D computes contour of quadratic function in 2D. % [X1,Y1,X2,Y2]=quad2d(A,B,C,win,interp) % % QUAD2D computes points on countour of 2D quadratic function % x'*A*x + B'*x + C = 0, for x in win. % % The computed curve(s) are returned in the format useful % for Matlab 'plot' function. % % Input: % A [2x2], B[2x1], C[1x1] parameters of 2D quadratic function. % win [1x4] defines function domains. % interp [1x1] number of lines used for interpolation. % % Output: % X1 [1x(interp+1)], Y1[1x(interp+1)] points on the 1st countour. % X2 [1x(interp+1)], Y2[1x(interp+1)] points on the 2nd countour. % % Example: % % A=-eye(2,2); B=[2 2]; C=[1]; % [X1,Y1,X2,Y2]=quad2d(A,B,C,[-1 1 -1 1]); % figure; hold on; % plot(X1,Y1); % plot(X2,Y2); % % See also L2Q2D, PQUAD2D, QTRANSF. % % Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac % (c) Czech Technical University Prague, http://cmp.felk.cvut.cz % Written Vojtech Franc (diploma thesis) 19.11.1999, 23.12.1999, 5.4.2000 % Modifications % 26-June-2001, V.Franc, comments repared. % 24. 6.00 V. Hlavac, comments polished. % default setting if nargin < 4, win=[0 1 0 1]; end if nargin < 5, interp=75; end % x-limits xmax=win(2); xmin=win(1); % y-limits ymax=win(4); ymin=win(3); %%%%%%%%%%%%%%%%%%%%%%%%5 a=A(1,1); b=B(1); c=A(1,2)+A(2,1); d=B(2); e=A(2,2); f=C; alpha=c^2-4*e*a; beta1=2*c*d-4*e*b; beta2=2*b*c-4*a*d; gama1=d^2-4*e*f; gama2=b^2-4*a*f; Ddx=beta1^2-4*alpha*gama1; Ddy=beta2^2-4*alpha*gama2; if Ddx >0 & Ddy > 0, yy1=(-beta2+sqrt(Ddy))/(2*alpha); yy2=(-beta2-sqrt(Ddy))/(2*alpha); if alpha < 0, ymin=max(ymin,yy1); ymax=min(ymax,yy2); else ymin2=max(ymin,min(yy2,ymax)); ymax2=min(ymax,max(yy1,ymin)); end end %%%%%%%%%%%%%%%%%%%%%%%%5 % curvature X1=[]; X2=[]; Y1=[]; Y2=[]; % go trough y value if Ddx < 0, % Ddx > 0 -> x=(-inf,inf) step=(xmax-xmin)/interp; for x=xmin:step:xmax, ay=e; by=c*x+d; cy=a*x^2+b*x+f; D=by^2-4*ay*cy; if D > 0, X1=[X1,x]; X2=[X2,x]; Y1=[Y1,(-by-sqrt(D))/(2*ay)]; Y2=[Y2,(-by+sqrt(D))/(2*ay)]; elseif D==0, X1=[X1,x]; X2=[X2,x]; if ay~=0, Y1=[Y1,-by/(2*ay)]; Y2=[Y2,-by/(2*ay)]; else Y1=[Y1,-cy/by]; Y2=[Y2,-cy/by]; end end end elseif Ddy < 0 | alpha < 0, % Ddy < 0 -> y=(-inf,inf) % or alpha < 0 & Ddy > 0 & Ddx > 0 -> ellipse step=(ymax-ymin)/interp; for y=ymin:step:ymax, ax=a; bx=b+c*y; cx=e*y^2+d*y+f; D=bx^2-4*ax*cx; if D > 0, Y1=[Y1,y]; Y2=[Y2,y]; X1=[X1,(-bx-sqrt(D))/(2*ax)]; X2=[X2,(-bx+sqrt(D))/(2*ax)]; elseif D==0, Y1=[Y1,y]; Y2=[Y2,y]; if ax~=0, X1=[X1,-bx/(2*ax)]; X2=[X2,-bx/(2*ax)]; else X1=[X1,-cx/bx]; X2=[X2,-cx/bx]; end end end % connect curvature if Ddx > 0 & Ddy > 0 & size(X2,2) > 0, if ymin > win(3), X1=[X2(1),X1]; Y1=[Y2(1),Y1]; end if ymax < win(4), X1=[X1,X2(size(X2,2))]; Y1=[Y1,Y2(size(Y2,2))]; end end elseif alpha > 0, step=(ymin2-ymin)/interp; for y=ymin:step:ymin2, upc=1; ax=a; bx=b+c*y; cx=e*y^2+d*y+f; D=bx^2-4*ax*cx; if D > 0, Y1=[Y1,y]; Y2=[Y2,y]; X1=[X1,(-bx-sqrt(D))/(2*ax)]; X2=[X2,(-bx+sqrt(D))/(2*ax)]; elseif D==0, Y1=[Y1,y]; Y2=[Y2,y]; if ax~=0, X1=[X1,-bx/(2*ax)]; X2=[X2,-bx/(2*ax)]; else X1=[X1,-cx/bx]; X2=[X2,-cx/bx]; end end end % connect curvature uindex=size(X2,2); if step > 0 & ymin2 < ymax, X1=[X1,X2(uindex)]; Y1=[Y1,Y2(uindex)]; end step=(ymax-ymax2)/interp; first=1; for y=ymax2:step:ymax, ax=a; bx=b+c*y; cx=e*y^2+d*y+f; D=bx^2-4*ax*cx; if D > 0, if first==1 & ymin < ymax2, first=0; Y1=[Y1,y]; X1=[X1,(-bx+sqrt(D))/(2*ax)]; end Y1=[Y1,y]; Y2=[Y2,y]; X1=[X1,(-bx-sqrt(D))/(2*ax)]; X2=[X2,(-bx+sqrt(D))/(2*ax)]; elseif D==0, if first==1 & ymin < ymax2, first=0; Y1=[Y1,y]; X1=[X1,-bx/(2*ax)]; end Y1=[Y1,y]; Y2=[Y2,y]; if ax~=0, X1=[X1,-bx/(2*ax)]; X2=[X2,-bx/(2*ax)]; else X1=[X1,-cx/bx]; X2=[X2,-cx/bx]; end end end end return;